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1. Introduction 

Gaussian wave functions have found increasingly wide application in recent years for molecular 
energy calculations [l]. 1 This is due to the relative ease of computation of the three- and four- 
center coulomb integrals; a circumstance that was first noted by Boys [2] and McWeeny [3]. Boys 
determined explicit formulas for all the relevant molecular integrals involving Is Gaussian functions 
centered at arbitrary points in space. He also observed that all other integrals could be obtained by 
differentiation of the resulting formulas. Harris noted that this procedure would be very tedious 
and he extended the explicit algebraic representation of the integrals to cases where the principal 
quantum number of the Gaussian function equals the azimuthal quantum number [4]. However, 
when the principal quantum number exceeded the azimuthal one the integrals were left as functions 
of differential operators. In this paper explicit formulas will be obtained for all cases by a different 
procedure than that used by Harris. 

Both the one- and two-electron coulomb integrals will be considered. The nuclear attraction 
integral is 

(Anlm | r- c ! | Bn'im') = j dr^^^Jr)^ (1) 

and the electron repulsion integral is 

(AnlmBn'l'm' | Cn"l"m!'Dn'"l'"m"') 

*.Hi^i)U^to*« ( 2<(f>ir', (2) 



■/■ 



where ri2 =| r ,— r2 | and r c =\r— C\. 

The usual expansion of the function r^ -1 in spherical [5] or cylindrical coordinates [6] is rejected in 
favor of the relation [7] 



277 2 J 



dk k~ 2 exp [ik • r], (3) 



This representation is favorable for the evaluation of coulomb integrals with basis functions of the 
following form 

<Wm£) = r A 2n + l Pi |m| (cos A )e im *A- ar A 2 (4) 

where (r A , A , $a) are spherical polar coordinates with origin at the point A and P; |m| (:c) is the 
associated Legendre spherical harmonic [8]. It will be seen that the use of eq (3) essentially 
reduces the task of evaluating the coulomb integral to the determination of the three-dimensional 
Fourier transform of ty n f m . 



1 Figures in brackets indicate the literature references at the end of this paper. 
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It is convenient first to reduce the dependence of the coulomb integral on molecular geometry 
by referring the coordinate axes at the atom centers to a new system with a common orientation. 
This was done by Harris. Such a transformation will be discussed briefly in the second section of 
this paper. By these means we isolate the basic integrals, whose evaluation is the main concern of 
this paper. 

The Fourier transform of tff n f m is readily found in section 3. The remaining integrations over 
the k variables are then completed for the nuclear attraction integral in section 3 and for the 
coulomb repulsion integral in section 4. The use of the Gaussian functions is briefly discussed in 
section 5. 



2. Transformation of the Integrals 

It was shown by Boys that it is convenient to expand the product tynin$n'i'm' m e( I s W an< ^ (?) at a 
point P on the directed line between A and B. Similarly, the product ifi^^ is expanded 
about Q. The polar axes of the coordinate systems at points P and Q are initially along B— A and 
D—C, respectively. The coordinate system at P is obtained from the one at A by rotation through 
the Euler angles (a, /3, y), collectively characterized as (PA), and translation without rotation a 
distance \P-A\=bRil(a + b). 

The basis functions in eq (2) are then transformed by means of the following expressions. 

The rotation of the spherical harmonic through the Euler angles a, /3, y is given by [9]. 

i 

Pi'"' (cos d)l im * = ^ Di mm \a, 0, y)P e lm ' 1 (cos d')e im '*'. (5) 

m'=-l 

When the spherical harmonic is translated without rotation a distance t along the polar axis, it is 
transformed [10] 

i /l+\m\\ 

r7V m| (cos 6)= Y tt-ir'JPjW (cos 0'). (6) 

Also, the trinomial expansion will yield 

r ^ 2n== f 4 + ^rRirp cos Op + y-^Rt 2 
\ a + b (a + bjr 

(7) 

-2 "\ j (-T^i)" 2 (^^T'rp^irp cos e P )»2. 
*-* n\\ni\n^\ \a-rb / \a-rb J 

(2»/=n) 

The use of eqs (5), (6), and (7) permit the expansion of the product ^nim^n'Vm' in terms °f tne 
coordinate system centered at P. Included in this expansion is a product of two Legendre poly- 
nomials which may be expanded [11] 

Pj\m\( x )p\m\( x ) = y C mm>pm>-m( x y (g) 

We also note that by multiple use of the recurrence formula [12] for the Legendre polynomials we 
can obtain 

x^tKx) = ±a n P^ n (x) (9) 

n=-j 

where a n = if n=— j + 2r+l for integer r. The a n can easily be determined for any particular 
case. A similar series of transformations is applied in the coulomb repulsion integral to transform 
the product ^JfT^^^^A the coordinate system is centered at Q. 

The last step in transforming the integrands for both the nuclear attraction and coulomb 
repulsion integrals is to rotate the coordinate axis of P in the first instance to have its polar axis 
directed along C— P and for the coulomb repulsion integral the axes at P and Q are rotated to have a 
common orientation with the polar axis Q—P. In both these cases the Euler angles for the rotation 
at P will be denoted by (P'P). These results and the definitions in the appendix allow us to factor 
out the dependence on molecular geometry from eqs (1) and (2). 
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3. Nuclear Attraction Integral 



The nuclear attraction integral can be written 



(4nlm\r c - 1 \BnTm')=e-*i J £ Di m "*(PA)Dv m '°'(PB) 



7 =-l a ' = -V 
I V 



j==\a\ j'=\cr'\ (In^n) (ln' = n') 
i i 1 



i l 

^sb^ & *r (-& «■) "t^ *r fe *r •"» 

ycji,-^' s aA - 2 Dj+^'-^iP'P) 

j k=-(n 2 +n' 2 ) M=-(J+k) 

(dTr p W+* n i+ n i )+ W n ' 2 Pj+k m (cos (^e^p-^+^pVeT 1 . 
It is readily seen that we must evaluate the integral 

/m = (drrp 2N+I PiJ^ (cos <bp) e**+p-*»irc- v (11) 

NL J 

The primes have been dropped in eq (4) but the coordinate system has its polar axis along C — P and 
its origin at P. Although M = in this coordinate system from symmetry considerations, we shall 
nonetheless retain the factor Pi) M \ (cos <j))e iM(i) since the more general expression is required for 
the electron repulsion integral. Substituting (3) into eq (11) we obtain 



NL J 

OU2 J 



2IP J ( 12 ) 

\dTpr P 2N + L Pi M (cos (f) P )e-^ r f> +i ^^H>- 

Without loss of generality <f)jc can be set equal to zero, so that 

exp [iA:-rp] = exp [ikrp (cos Op cos i//-!-sin dp sin i// cos (/)/>)] 

where i// is the angle between the CP axis and k. Consider 

J=\ dr P r P 2N+L+2 \ \ dS P d<t>p sin 6 P P L M (cos ftp) 
Jo Jo J 

exp [iA;rp (cos 0p cos t// + sin 0p sin i// cos </>p) + iMcfrp — U 2 r P 2 ] (13) 

and integrate over the angular coordinates to obtain 

J = {2tt?I 2 j^ 2 P L M (cos *) rdrprp^^VL+i^Me-^. 

Utilizing the Weber and Sonine formula [13] 



/; ^(a,)e-->-.=^A^ lfl (ft±£ ;„+!;-£) (14) 
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we obtain finally 

J=J ¥ iL nL+tmv^ 1 ^ P,M (cos * *^ N+L + 3 / 2 ; *+3/2;-J^). 

Substituting 7 into eq (12) and setting M = O y we obtain 

. ^ r(7V + L + 3/2) 

7 ^ 2 l+1 tt 1 / 2 T(L + 3/2)f/ 2Ar+2/ ^ 

J dkk L ~ 2 P L (cos </0iF, (iV+Z, 4- 3/2; L + 3/2; - ^) exp [i£ • CP]. 
Integration over the angles yields 

7 " L = jS\^\ dkkL ~ ll2 J^l^ k \ CP \)^ (#+L + 3/2;L + 3/2;-^ 



with 



C = 2- f ^_ IW+I + 3/2) 

T(L-f 3/2)i/ 2A7+2L+3 



This integral, considered as a Hankel transform [14], is given as 



F„, = 1^ /°21 T]2\rp\2 

\CP\ L+1 U 2N+2L+3 23 



1, L+3/2 
L+l/2, JV+L+3/2, 0, 



where the G function is a polynomial in generalized hypergeometric functions [15]. 

An alternative formula may be obtained by relating the confluent hypergeometric function 
to the Laguerre polynomial [16], 

lfl (^ + L + 3/2;L + 3/2;-^)=^-e-/-^ 1 / 2 (|,) 

and(a) L = (a)(a+7). . . (a + L-1). 

The Laguerre function is given as the power series [16] 



L ^ I/2 te/~§ HS l N . s ](27/FS! 



and therefore 



C N\ A /,v+i+i/2\ i 

\dk k^-^jL+miklCPDe-W 4 " 2 (15) 

Jo 



JV! _ \CP\ L A Hf ^+t+i/2 
IXL + 3/2) !/*»+« ^# S! 



s=o 



r(L+S + l/2)exp(-l^|CP| 8 )i J Pi(-S+l;L+3/2;£/ 2 |C/ > | i! ). 
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4. Coulomb Repulsion Integral 

The coulomb integral, in the most general case, will reduce to 
[AnlmBn'l'm'\Cn"l"m"Dn , "l"'m'"] = e- x i- x 2 



l V l" l" 



X 2 2 2 Dr T *(PA)Df°XPB)Df/°"*{QC)Df"°"XQD) 



I V I" I' 



v v v v ( bRl Y" j ( - aR i Y~ r ( gggg V"" j " ( -CRi T'-y 
i+\a\\/r + \<T'\\(i"+\<T"\\(r" + \<T'"\\ 

j+\<r \ \j"+ <r'A/'+-o-"/\7"" + W"\) n ** „^ «„ £ „,„ P „, 

j i |/ \j \ \/ \j i |/ w i |/ (I^= n ) (2n.=n') (Xn'/=w") (ln'."=n'" 



i i 



n\n'\n"\n'"\ /26gA"2 /_2oRi\»i 

» 1 !n 4 !7i 3 !ni!^!nyre'{!n^!n'i!ii"{!n"^!re"^!Va+6/ I a + b) 

/2dR 2 \"'i ( 2cR 2 \"«" / 6/?. \ 2 "3 / a«,\ 2 "3 / <#? 2 \ 2 "s / c«, \*"a" 



\C+d) \ c + d) \a + b) \a + b \c + d \c + d 



n 9 +n' n"+n'" 

2 2 2 (J+k)- 



./ J' A:=-(n 2 +n£) A-' = -(n£'+nf>") M=-(7+*)- 

J+k \cos v Pl )rj, +k ,{cos u Q2 )e vv q 2 ^p\' p\ Q2 

ri2 

where (J+k)- is the smallest of (J + k) and (J' + k'). Evidently we need to evaluate the follow- 
ing integral 

I M nln>l> = \ dT P1 dTwr P ™+LP L M(cos M V v ' + / ' V'(cos 6 Q2 ) 

eiM(<t> Q2 -<t> P1 )-u2r 2 pi -V*r 2 Q2 (17) 

where the polar axis is along Q — P. Again using the transformation of eq (1) we can separate 
the two-electron integral into a product of integrals each over only the coordinates of one electron. 
These integrals are, of course, very similar to the nuclear attraction integral that has just been 
evaluated. Therefore, we can immediately write down the integral over k. 

[m nlnili = c [ dk k L + L '- 2 P L M (cosilj)P L !"(cos ifj)eWP-Q) 

j (18) 

.F^N + L + Z^L + m^-^IW^^N'+L' + m^L' + m'.-k'IW 1 ) 



with 



C = i _ )L . 2 -L-V- llri L + V. 



[/2iv + 2L + ayw + tu + 3Y(L + 3/2)r(L' + 3/2) 
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The angular integrations are readily performed using eq (8) and the relation [17] 



e ikrcos* = ^ (2/ + 1) i l ji(kr)Pi( cos i//)- 

1 = 



(19) 



Equation (18) reduces to 

I M NLN .,, = ^C^VC^ JLL , f" dk W+vjAk\PQ\) 

J ■> 

iFi(N+L + 3/2; L + 3/2; - *»/4DVi(JV + U + 3/2; L' + 3/2; - F/4F 2 ) 

J 

The confluent hypergeometric functions can again be given in terms of the associated Laguerre 
functions. We then find 

IU) = ^K} f § ( _ )S + s.(K+L+i\/N'+L> + h 

(J! (L + 3I2ML' + 3I2) N .£ Z^ > \ N - S J\ N'-S' 

2-us 



->-US + S') f « 1-2/ 172 _|_ T/2\ 



U 2 + J 72 
For y 2 = y-T 2T72 , we again utilize the Weber and Sonine formula to obtain 

r . n v AW! 2^'~ l \PQ\ J 

1(J) = IT* 



(L + 3/2 UL ' + 3/2 V T(J + 3/2) 

y y ^r^^^Ym,, 



S=0 S'=0 



r 

y- 25 ^-"' „ / L + L'+7 + 2(S + S')+I |W\ 
SIS'! lFl I 2 ;J + W; —)• (20) 

The substitution of (20) into eq (19) completes the derivation of an explicit algebraic formula for the 
coulomb repulsion integral. 

5. Discussion 

The connection with the usual expression of these integrals in terms of error functions and 
derivatives of error functions may be seen from the representation of the error function by the 
confluent hypergeometric function [18] 






There is little reason, however, to replace the confluent hypergeometric functions with the more 
familiar error functions. For computational purposes the confluent hypergeometric functions are 
quite adequate. 
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These formulas make it possible to extend the use of Gaussian basis functions to higher 
principal and azimuthal quantum numbers which is required if the Hartree-Fock limit is to be 
approached for most molecules. The integral computation time will not materially increase as the 
quantum number is increased; the major difficulty is the inherent inaccuracy of the Gaussians 
(i.e., the wrong cusp values and asymptotic behavior) and a large number of basis functions is 
required to overcome this deficiency. For example, for a simple system like H2 about twice as 
many Gaussian basis functions as Slater functions are required to achieve comparable results [19]. 
Therefore, the main problem remains the handling of the very large blocks of numbers that are 
generated in the Hartree-Fock procedure. 

6. Appendix I. Definition of Terms 

aA + bB 



Ri 


= \A-B 


Rz 


= \C-D 


X, 


abRi* 


a + b 


x 2 


_cdR 2 2 



p= 



Q = 



a + b 

cC + dD 
c + d 



c + d U 2 = a + b 

V* = c + d 
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